Geometric structure of coexisting phases found in the 
Biham-Middleton-Levine traffic model 
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I. INTRODUCTION 



Raissa M. D'SouzeQ 
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The Biham-Middleton-Levine traffic model is perhaps the simplest system exhibiting phase tran- 
ce , sitions and self-organization. Moreover, it is an underpinning to extensive modern studies of traffic 
q_) . flow. The general belief is that the system exhibits a sharp phase transition from freely flowing 
fj^ ' to fully jammed, as a function of initial density of cars. However we discover intermediate stable 
phases, where jams and freely flowing traffic coexist. The geometric structure of such phases is 
highly regular, with bands of free flowing traffic intersecting at jammed wavefronts that propagate 
smoothly through the space. Instead of a phase transition as a function of density, we see bifurcation 
points, where intermediate phases begin coexisting with the more conventionally known phases. We 
show that the regular geometric structure is in part a consequence of the finite size and aspect ratio 
of the underlying lattice, and that for certain aspect ratios the asymptotic intermediate phase is 
on a short periodic limit cycle (the exact microscopic configuration recurs each r timesteps, where 
t is small compared to the system size). Aside from describing these intermediate states, which 
previously were overlooked, we derive simple equations to describe the geometric constraints, and 
predict their asymptotic velocities. 
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Modeling vehicular and Internet traffic, thereby gaining an understanding of congestion patterns and jamming 
phenomena, is an extremely relevant problem, with obvious practical ramifications. One popular approach is the 
I ' use of simple discrete cellular automata (CA) models, which capture aspects of the dynamics of discrete vehicles or 
, packets. One of the most cited examples of such a CA model is the Biham, Middleton, and Levine model (BML) of 
t-h ' two-dimensional traffic flowQ. At the time of this writing, Ref. has received over two-hundred citations in the 
| scientific literature, and it serves as a theoretical underpinning for the physicists' approach to modeling traffic. Note 
that using techniques from physics to model traffic has been a fruitful research area for more than a decade, and 
continues to be. For recent reviews see Refs. QEiOLS- 

The BML model describes two species of "cars" moving on a two-dimensional square lattice, with periodic boundary 
conditions. The model is extremely simple, yet the behaviors it displays, extraordinarily complex. The system shows 
what appears to be a phase transition from having all cars freely moving at all time steps, to complete gridlock, 
where no car can ever move again. In addition, in all the phases, the system becomes fully correlated, forming 
a range of interesting stable self-organized patterns. It is perhaps the simplest model where one can study both 
phase transitions and self organization. This model has recently become increasingly of interest to the combinatorial 
mathematics community, as it continues to elude rigorous theoretical analysis 0. 

We implement the BML model and experimentally study its behaviors. We discover stable intermediate states that 
have never been reported before for the BML model, with highly structured geometric patterns of wave fronts of 
jams moving through otherwise freely flowing traffic. We show the geometry of these patterns arise due to the finite 
size and periodic boundary conditions of the underlying lattice. We also show that the aspect ratio of the lattice 
imposes geometric constraints which restrict the patterns, and derive simple equations describing these geometric 
constraints which allow us to calculate the asymptotic velocities. For certain aspect ratios we can prove that the 
intermediate configurations end up on a short periodic limit cycle — the exact microscopic configurations recur each 
t timesteps (where r is small compared to the system size) — hence these states are stable for all time. For the other 
aspect ratios, we show the intermediate states are at least metastable, lasting as long as we could simulate them. By 
establishing the existence of these intermediate states, we show that the conventional beliefs about this model need to 
be reexamined. Contrary to the evidence published elsewhere, only on smaller spaces do we see evidence for a sharp 
transition from freely flowing to fully jammed configurations as a function of the initial density of cars, p. Instead we 
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FIG. 1: Typical configurations observed for the BML model on an L x L system of size L = 256. (a) The free flowing stage, 
where all particles advance during each update (v — 1). Note the ordered stripes of alternating east and north-bound cars. 
The width of the stripes increase, on average, with density p while in the low density phase, (b) A fully jammed configuration, 
consisting of one global jam. Note the jam length, V2L, is larger than the system size, (c) A high density, random jam 
configuration. 

observe bifurcations as a function of p, where different phases can begin to coexist. The bifurcation points, the range 
of the windows for phase coexistence, and the number of coexisting phases, depend on the size and the aspect ratio 
of the underlying lattice. Considering the amount of ongoing work on this model, and its use in large scale, complex 
simulations of traffic, we believe it is important to understand these new observations. 

This manuscript is organized as follows. In Sec. [H] we review the BML model and relevant past work. In Sec. IIIII 
we describe our simulations and empirical results. Section llVl contains a discussion of the kinetic pathways, geometric 
constraints, and derivation of the velocities for the intermediate states. Finally, in Sec. El we summarize and discuss 
open questions and areas for further inquiry. 

II. THE BML MODEL 

Consider two species of particles (i.e., "cars"), east-bound and north-bound (which we also interchangeably call 
"red" and "blue" respectively), which populate a two-dimensional square lattice with periodic boundary conditions. 
Each lattice site can be in one of three states: empty, occupied by an east-bound car, or occupied by a north-bound 
car. The cars are initially distributed uniformly at random along the lattice sites, with spatial density, p (usually 
taken to be the same for both north and east-bound cars). The discrete time dynamics has two phases. On even 
steps, all east-bound cars synchronously attempt to advance one lattice site towards the east. If the site eastward 
of a car is currently empty, it advances. Otherwise, it remains stationary (even if the eastward site is to become 
empty during the current time step). On odd time steps, the north-bound cars follow the analogous dynamics, only 
attempting to advance to the north-bound site. The dynamics is fully deterministic. The only randomness is in the 
initial condition. Furthermore, the dynamics conserves cars, and does not allow for an east-bound car to change its 
row, nor for a north-bound car to change its column. So on a L x L lattice, there are 2L conservation laws. 

If initialized with a low enough density of cars, the system eventually self-organizes into a configuration where all 
cars can move at each time step (each car has asymptotic velocity equal to unity). A typical such configuration is 
shown in Fig.^a). If initialized at slightly higher density, the cars are blocked by other cars, until eventually all cars 
end up participating in one large global jam, where no car can move (asymptotic velocity equal to zero). A typical 
global jam is shown in Fig. ^b). The transition between the two behaviors has been thought of as sharp, showing 
characteristics of a first order phase transition. Initialized at much higher densities, small jams begin simultaneously 
throughout the lattice and merge almost immediately with other small jams, leaving all cars blocked (with all velocities 
equal to zero). In this high density phase, the system has no time to self organize, and instead of one global jam, we 
observe a collection of small random jams. This latter type of jam has been compared to traffic in a large city during 
"rush hour": a car might escape one jam, only to quickly join the tail of the next. We expect that for an infinite size 
system, the fully jammed state resembles this random type of jam. 

Most of the understanding of the BML model has come from numerical simulations^, 0, IE El El El- 
Theoretical analysis has been limited to mean- field approaches l?l IT3. Il5l Il6| . and attempts to start with continuum 
hydrodynamic equations and formulate an equivalent discrete model|l7|. There are general beliefs about this model, 
that the transition is first order and that the critical density, p c , decreases with increasing system size, possibly 
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reaching the value p c = as the system size approaches infinite |l(. The BML model has been simulated extensively, 
but there are inconsistencies in the literature and lack of detail of numerical implementations (such as the size of the 
ensemble beingaveraged over). Details of numerical studies have been published only for small systems, on the order 
of L = 10 — 50[7, 8] . Larger systems have only been studied coarsely, or in the context of self-organized versus random 
jams@,03. 

Despite extensive numerical simulation, the existence of stable intermediate phases (with < v < 1) has not been 
explicitly reported previously. Fukui and Ishibashi[8j do show evidence of the existence of an intermediate phase in 
one plot. They note that for intermediate values of p, v "fluctuates around a certain value for a long time" . The value 
shown in their plot is extremely close to the values we observe (plotted in Sec. Aside from this comment, they 
do not pursue the issue further. A careful study by Torok and KerteszfTlj contains precise details of their numerical 
simulations. They are studying a variant of BML with faster convergence times (called the green wave model). Since it 
is not possible to theoretically predict the convergence time, they estimate it, and apply the following, very reasonable, 
empirical heuristic. If a realization has not reached a state, with v rs or v rs 1, within an allotted time (taken to be 
5 times the estimated convergence time), that realization is discarded. We can only assume some of the studies of the 
BML model may have used a similar criteria of discarding "non-converged" states. Note that for continuum models, 
intermediate phases of jammed wavefronts moving through otherwise freely flowing traffic have been reported. See 
for instance Ref. [Is| . 

Sensitivity of the BML model to boundary conditions has been reported previously. Martinez, et a/.|l2| study the 
dilute limit (p — > 0), and show that different results are obtained for an "entangled" torus versus a conventional one. 
They raise interesting questions about how to get at the bulk properties using only finite size simulations, but do not 
quantify nor pursue the effects further. Chau, Wan, and Yan 13| study the BML model on the torus with random 
boundary conditions (BCs) (meaning, particles moving off the right (top) edge reappear at some randomly selected 
site on the left (bottom) edge). They claim v > whenever p < 1, and hence dismiss such systems as being "not very 
interesting" . They also note that the velocity and critical density depend sensitively on the choice of BCs, but they 
also do not pursue the effects further. 

III. SIMULATION RESULTS 

We implement the BML model on square lattices of finite size L x L', for a range of sizes and varying aspect 
ratios. For square aspect ratios (i.e., L — L'), the lengths studied range from L = 64 to L = 512. We also study 
rectangular aspect ratios where the width, L, is an integer multiple of the height, L' , and where L and L' are relatively 
prime. On each lattice, we implement a range of densities p, studying at least ten realizations for each density. All 
simulations are implemented on a special purpose cellular automata machine, the CAM8 19] . The CAM8 performs 
approximately one-billion site updates per second, comparable to a modern high-performance desktop computer. The 
main advantage of the CAM8 is that it allows for excellent visualization of the system, with no overhead in the 
rendering, giving us live video output of the dynamics of the system. As discussed in Sec. IIVI visualizing the kinetic 
pathway gives crucial insight into the formation of the intermediate states. 

The main subtlety involved with simulating this model is determining the convergence time (i.e., the time is takes 
to reach the asymptotic behavior). All realizations were simulated until converged (v = 0, v = 1, or the periodic limit 
cycle was reached) or for times out to at least t — 2 x 10 6 time steps. Of the realizations that had not yet converged, 
many were simulated for orders of magnitude beyond. We find this a reasonable compromise, since the compute power 
to simulate all samples to times greater than 10 8 is beyond our current capacity. Throughout the remainder of the 
manuscript, we refer interchangeably to the east-bound cars as "red" and the north-bound as "blue". 

A. Square aspect ratios 

For small size systems we actually observe the predicted behavior of a sharp transition from freely flowing to total 
gridlock. Figure (0}) is for an L x L system with L = 64. It is a plot of the final average velocity observed for each 
realization, v, versus the density for that particular realization. Note we are plotting the average velocity for each 
individual realization, not an average over all realizations initialized with the same p (hence error bars are on the 
order of the size of the plotting symbol used). Surprisingly, when we implement systems with L > 64, we observe a 
bifurcation where two phases start to coexist, as we go from low to intermediate values of p. The second phase that 
emerges has average velocity v ~ 2/3, as shown in Fig.©. In addition, these intermediate states have a very well 
defined geometry, of bands of red stripes with slope one-half, criss-crossing bands of blue stripes with slope two. An 
example of the geometry is shown in Fig.Q. Jammed wavefronts are located at the intersections of the bands, and, 
as the systems evolves in time, move as solid structures uniformly down towards the southwest with unit velocity. 
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FIG. 2: Plot of the average velocity for each individual realization, v, versus the density, p, for that realization, for an L x L 
lattice: (a) L — 64, (b) L = 128, (c) L = 256, (d) L — 512. Note the emergence of an intermediate phase for L > 64. The value 
of v for the intermediate state becomes more crisply denned with increasing system size, and that the window of coexistence 
between the phases broadens. For the L — 512 system, the average value of v in the intermediate phase is (v a ) = 0.673 ±0.005. 
The dotted line shown in (b)-(d) is the prediction from Eq.l(24j. 



Particles are freed from the head of each jam, but a like number of new jammed particles aggregate at the tail. As 
discussed in detail below, the underlying lattice imposes constraints on the allowed topologies of the configurations. 
In Sec. IIVI we derive a simple formulation of the constraints. 

All realizations contributing to Fig.(J2J) were simulated to at least ( = 2x 10 6 full updates of the entire space, with 
various realizations simulated for t > 10 8 . We do occasionally observe a realization persist in the v ~ 2/3 state for 
orders of magnitude, then suddenly jump to either v — or v — 1, with the latter being more common. Regardless, the 
intermediate state is at least metastable, persisting for longer than we could simulate most realizations. Furthermore 
it is "universal" in the sense that the value of v ~ 2/3 is independent of system size and density, and all systems that 
do not go to v = or v = 1 go to the same intermediate state (i.e., same geometric structure and approximate value 
of v). The reason we do not average over the individual realizations, is that it would obscure the behavior. Instead 
of displaying the three distinct quantized states, averaging would produce a deceitfully smoothly decaying curve. 

From our data, it is difficult to determine the exact bifurcation point where the phases begin to coexist, and the 
point where they cease to. We attempted to identify the factors that distinguished realizations which converge to v = 
from those, with the same density, that converge to v ~ 2/3. We first investigated connections to anisotropy, such as 
an imbalance between the total number of red versus blue particles. But we found no correlation. The probability 
a realization would jam or go to the intermediate state is independent of this asymmetry. We also looked at a more 
fine-grained measure: the line density of red particles versus the blue. Again we found no correlation between this 
asymmetry and the likelihood of jamming. 

In Sec. II VI we discuss the kinetic mechanism observed, which gives rise to the interleaved band structure exhibited 
by the intermediate states. As mentioned above, the jam interface moves ballistically, with unit speed, towards the 
southwest. The width of the jam interface can fluctuate. It seems for the L = 64 system, the fluctuations are large 
enough that eventually the head of one jam meets the tail of the previous, continuing until eventually one global 
jam forms. Kinetic mechanisms aside, it could be that the convergence times for L = 64 are short enough that the 
intermediate state cannot be considered metastable. Fig.Q shows the median convergence times observed for the 
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FIG. 3: Examples of typical intermediate geometry, (a) A system with square aspect ratio where L — 512, ui r = 1/2 and 
ujb = 2. (b) A system with rectangular aspect ratio where L — 377, L' = 233, uj,- = 1, and uib = 3. Note that in (a) there are 
many disordered, random cars in the space between the bands, yet in (b) all cars are ordered. We find this crisp order shown 
by the latter example for all realizations studied on rectangular aspect ratios with L and L' relatively prime. If L and L' are 
not relatively prime, random disordered cars located between the bands persist. Note that Fig.JHJ is a close up of the region 
that has just shed from the jams in (b). 
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FIG. 4: Median convergence times for the different square aspect ratios simulated. These results are consistent with previous 
studies. Note a value of t = 2 x 10 6 should be interpreted as t > 2 x 10 6 . 

systems with square aspect ratios. Note that the a value of t = 2 x 10 6 really means t > 2 x 10 6 . 

B. Rectangular aspect ratios 



We also implement the BML model on systems with varying rectangular aspect ratios. In particular we study 
lattices where the two lattice lengths are subsequent Fibonacci numbers. Figure (JSJi) is the plot analogous to those 
shown in Fig.JSJ), for a system with L — 89 and L' — 55. Note we see the same intermediate velocity of v ~ 2/3, but 
we also see the emergence of one more possible phase with v ~ 2/5. All the configurations with v ~ 2/3 resemble the 
one shown in Fig. 103), with one red band wrapping around the f-axis, and three blue bands. The jam points are at 
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FIG. 5: (a) The average velocity for each individual realization, v, versus the density, p, for that realization, for anixL' lattice, 
with L = 89 and L' = 55. Note the appearance of yet another well defined intermediate state. Note also, the bifurcation point, 
where phase-coexistence ceases, has not yet been reached, despite the range of p > 0.44. The dotted lines are the predicted 
velocities from Eqs. ll8H and 11211 . v r i ~ 0.7430 and v r 2 ~ 0.3707. The empirically determined average values are respectively, 
(v r i) = 0.700 ± 0.002 and (1V2} = 0.364 ± 0.004. The geometries of the two types of intermediate states are distinct: (b) A 
typical configuration with v ~ 2/3. Interface slope s — 1; (c) A typical configuration with v ~ 2/5. Interface slope s w 2/3. 



the intersections of the bands. Note the crisp, regular geometry. Figure is for an equivalent, but larger, system 
with the same Fibonnaci aspect ratio but with L = 377 and L' = 233. This figure more clearly illustrates the highly 
ordered geometry. It also includes the definitions of several of the parameters used in the analysis in the subsequent 
sections. All realizations with v ~ 2/5 resemble the one shown in Fig-JS^), with one red band and approximately two 
blue bands (though the latter are not so clearly defined). 

One of the most striking differences when comparing these rectangular aspect ratios to the square, is that for the 
rectangular, the intermediate configurations are exactly periodic with a short recurrence time: the exact microscopic 
configuration of particles repeats every r updates. We observe systems of size L = 89 and L' — 55, settling into 
the periodic behavior typically in a time less than t ~ 100, 000 updates, with a period on the order of r ~ 6000 
updates. Figure is a plot illustrating such typical behavior. Another striking difference, is that not one realization 
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FIG. 6: For simulations with L = 89 and L' = 55 and p = 0.38, we plot the sample number versus, (a) the time to reach the 
periodic limit cycle; (b) the period of the cycle, r. Eight of the ten realizations simulated reach final configurations of type I, 
as shown in Fig.lJ^J)). They all have the shorter values of r. Surprisingly six of these eight realizations have the same period, 
t = 5114, though their microscopic configurations differ. The other two realizations, with significantly larger values of r, are 
of type II, as shown in Fig.J^t). 



7 































o ! 




X 


X 


XX | 


X 


X 




X 


X 


X 


o] \i) 


X 


X 




X | X 




X 


X 


X o | 


X 


X 


x ;o ; 


X 


X 




X 


X 


X 





_J?l_i_. 






X 


°j°Li„L, 


X 


X 




X 


X 





o ! ' ' 


X 




X 






X 





°J °J J 1 


X 





_°J?J_i_i_. 


X 








! ! ' 


X 





































1 1 1 











' * 1 






t=0 






t=l 






t=2 






t= 


3 


























i i i 



























| X 


X 





|x | [x 


X 




X 


Pal 1* i I L " l X I 


|_x 




X 




X 









'X 







X iii 




X 


p 






X 





1 







X 





iO 








i : i i 










t t r r 






1 


=4 








1=5 



X o 

X 



o 



o ! ' ! 



X o 

1 x^ 



FIG. 7: Shedding from a jam. 



has jammed fully (v = 0), despite the fact that the largest density simulated is p = 0.45. Note that for the square 
aspect ratios, the bifurcation point where the phases cease to coexist is at p ss 0.40. Another striking difference is the 
lack of disorder for the relatively prime systems with rectangular aspect ratios. In Fig.(P^,) there are isolated particles 
("dislocations") moving in the area between the bands. However for the rectangular case, all the particles manage to 
join the ordered bands. 



IV. GEOMETRIC CONSIDERATIONS 



A. Kinetic pathways 



By watching the dynamical evolution of the system, starting from the initial configuration, the mechanism by which 
the intermediate phases form can be observed. Often one global jam initially begins to form, yet the head of the jam 
just fails to meet up with the tail, leaving a few lattice sites of distance between them. Particles shed from the head 
as soon as allowed by the local environment (since all particles move whenever possible), leaving with a well defined 
order. See Ref. [2£j], for video images of the dynamical evolution. 

To understand the pattern formed by the shedding, consider first one row of a solid isolated block of red particles. 
Since the particles only advance provided the site they wish to occupy is empty, a particle would leave the head of 
the jam only every other timestep. However, we can now consider a diagonal interface formed by a triangular block 
of red particles in contact with a triangular block of blue particles. See Fig. 10. Note red (eastbound) particles are 
represented by "x" and blue (northbound) particles by "o" . Each subsequent time step illustrated corresponds to one 
complete update of the space (i.e., one north-step followed by one east-step). Recall, all particles of the same species 
update synchronously. Step 1: No "o's" are able to move. But, during the east-phase of the timestep, the first "x" 
moves away from the jam, opening up a space. Step 2: The first "o" advances. This blocks all other "x's" in the 
original row, yet opens up a space for an "x" , one site from the original "x" along the southwest diagonal, to advance. 
The original "x" to move, also continues advancing. Step 3: That first "o" continues advancing and opens up a space 
for an "x" in the original row to advance — two time steps delayed from the first "x" to move in that row. The "o" 
below this original "o" is currently blocked. However an "o" one site from the original along the southwest diagonal 
is now free to advance. Note this "o" blocks the "x" in its newly occupied row from advancing, yet opens up a space 
for the "x" one row southward to advance. Hence a pattern emerges: a red particle sheds from within the same row 
only every third timestep, yet from a site further southwest every timestep, yielding red bands of density p r = 1/3 
with slope s r — 1/2. Likewise for the blue particles, a blue particle sheds from within the same column only every 
third timestep, yet from a successive site along the southwest every timestep, yielding blue bands of density pb = 1/3 
with slope s;, = 2. The jams occur at the intersections of the bands, and the interface of a jam has slope s. Typically, 
s = 1. In Fig.Q one can see the order beginning to emerge. 

More illuminating is to view a closeup near the jammed regions. Fig.© is a zoomed in view of the region near one 
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FIG. 8: A zoomed in view of a jam shown in Fig.ll^b). 



of the jams shown in Fig. 103). Note the order that exists in the region above the jam interface, which has shed from 
the jam in the manner described above. In this region there are alternating diagonal stripes of blue, red, and empty. 
The stripes have "phase locked" — on the next update the blue stripes move into the empty stripe regions, leaving 
room for the red stripes to move on the subsequent update, and so on. Hence the system has organized itself into the 
highest density packing that still allows all particles to move with v = 1. 



B. Winding number 



Since the system lives on a torus, the bands must wrap seamlessly around it. Noting these facts, we can develop 
a mathematical expression for the number of red bands, io r and the number of blue, ujb, (referred to respectively 
as the "winding" number for red and blue), that must be present. Essentially, we can calculate the length of the 
regions of slope-1/2, slope-s, and slope-2 that must be present for the bands to wrap around the torus. (The slopes 
are respectively those of the red band, the interface, and the blue band) . Consider an L x L' lattice, and a red band 
starting in the lower left hand corner. Moving out along the ^-direction, there are k sites with slope-1/2 and k' sites 
with slope-s. Similarly for a blue band starting in the lower left hand corner, there are m sites with slope-2 and m' 
sites with slope-s. The constraints are: 



k + k 



— k + sk 

m + m 
2m + sm 
k 



L 

u r L' 
L 

u b L' 
m' , 



(1) 

(2) 

(3) 
(4) 
(5) 



where, due to the lattice, {k, k' ,m,m'} are positive integers, and {u> r ,u>b} are either positive integers or equal to l/g, 
where g is a positive integer (e.g., for the configuration shown in Fig.(|3]}), u> r = 1/2). The final equation, Eq. J5J, 
expresses that the length of the jam interface must be the same for the red and blue cars. Combining this system (of 
five equations in seven unknowns), we can solve some combination of variables in terms of the others. Solving first 
for k!\ 



k' = (2to r L' - L) /(2s - 1). 



(6) 



9 



L/L' 


s 






V 


1 


1 


1/2 


2 


0.673 ± 0.005 


1 


1 


1 


1 





2 


1 


1 


4 


~ 2/3 


2 


1 


2 


2 





5/3 


1 


1 


3 


~ 2/3 


3 


1 


2 


5 


~ 2/3 


(l+V5)/2 


1.17 


1 


3 


0.700 ± 0.002 


(l+VI)/2 


7/10 


1 


2 


0.364 ± 0.004 



TABLE I: Allowed winding numbers, u) r and ujb, for lattices with different aspect ratios, L/L' , and interface slope s. The values 
of ii which include error bars are from our numerical simulations. 




FIG. 9: A typical realization with L — 769 and L' — 256 (i.e.. L — 31/ + 1, where the additional lattice site is to make the 
lengths relatively prime). Note, that as predicted in Table HJ the configuration is consistent with s = 1, u> r = 2, and uib = 5. 
Realizations with L = 768 and L' — 256 (i.e., L — 3L') have a similar structure, but surprisingly lack the crisp order, and 
instead have disordered cars at random locations between the bands. For images of the latter, see Ref. |2(J. 



Using this, we can solve for LUf, in terms of oj r , s, L, and L' , and obtain: 



2L (2 -a) ( n L 
UJ b = T7- 7777 — TT l 2 "r - 77 • (7) 



U (2s - 1) V L 

Knowing that the interface slope s ~ 1, we can tabulate the allowed values of u r and lo^ in terms of the aspect ratio 
of the space, L/L' . The allowed values for various aspect ratios are listed in Tabled Note, the value u> r = 1/2 means 
the red band has only reached height L'/2 in traversing distance L (as in Fig.J^Ji)). We implement systems with the 
various aspect ratios shown in Table [IJ and find the experimental configurations observed all match the predicted 
behavior. See in addition to the previous figures, Fig.©- The final lines in Tableware for systems where we observe 
empirically the values of u) r and u>b, and using Eq.(|7J) can predict the value of s. Recall u> r and u>b must be integers. 
The system seems to tune the value of s to allows this. For instance, if L and L' are successive Fibonacci numbers 
(i.e., L/L' = (1 + \/5)/2), a rearrangement of Eq.(jZJ) predicts s = 1.17, which matches our empirical observations. 
Note in Fig.J^t), the upper jam has small glitches where s > 1. 



C. Average velocity 

Figure|31 illustrates the typical geometry for realizations with v ~ 2/3. In this example, uj r = 1 and Wf, = 3. We 
label on the figure the lengths a, b, c, and d, denoting respectively the width of the red band, the blue band, the blue 
jam and the red jam. We define two discontinuous functions which will simplify notation later on: 

e r K) = { 1 / ifw ^ 1 (8) 

l/w r if w r < 1, 
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In other words, © r is the number of independent red bands in one column of the lattice. Q b is the number of 
independent blue bands in one row of the lattice. If we denote the density of the red band as p a , we can express the 
average number of red particles in a column of the space: 

e r (oj r )p a a = pL'/2, (10) 

Likewise, the average number of blue particles in a row: 

e b (uj b )p b b = pL/2. (11) 

Recall p is the overall particle density (red plus blue). Empirically we determined that p a = p b = 1/3, which is a 
consequence of the dynamics described in Sec. IIV Al Using the basic equations described in Sec. II V Bl and knowledge 
of the "typical" geometry, we can solve for the velocity of the intermediate state. Unfortunately, we have to consider 
the square and rectangular aspect ratios independently, and do not have one equation that describes all cases. The 
assumptions described above are valid for the rectangular aspect ratios, but fail to capture the full behavior of the 
square ones. 



Rectangular aspect ratios, type I: 



Configurations such as the one shown in Fig. J2J) have a regular, highly structured, geometry that is well described 
by the formalism in Sec. IIV Bl We find that the values of k' predicted by Eq. exactly match those obtained 
experimentally. Typically there are rij jams, with the particles equally divided amongst them. The interface width 
per jam is k[ — fe'/njt, where n& is the total number of jams one red band is involved in, as it wraps once around the 
x-axis. See Fig.(|3Ji), for example, where ~ 2, and rij — 3. 

The structure of the jams and the relevant geometric factors are shown in Fig. HOf a). The jams form trapezoidal 
shapes of width 6. From simple geometric considerations, we can show that S = k[ sin(# — <p) = Kk[. The jam width 
per line, d, and per column, c is 



= d = 



sm 



cos£ 



cos£ 



(12) 



However, note this is only true provided that there are enough particles available in each column and in each row. 
Otherwise, 



d = rf max = pL/2uj r Q r (uj r ), and c = c max = pL'Q b (uj b )/2uj b . 
From similar geometric considerations we can show the length of the blue jam, l b , is approximately 



b I cos < 



tan((7 — 



Similarly, the length of the red jam, Z r , is approximately 



a I cos £ 



sin £ 



67. 



atp. 



(13) 



(14) 



tan (a — q 

The total number of particles involved in the jams, J, is the number of jams multiplied by the width and length: 

J sa rijS (l b + l r ) 



(15) 



(2uj r L' -L)k 



P L/e b (uj b ) 7 + -pL'/e r (u r ) <p 



J 



(16) 



For the data plotted in Fig. EJ L' and L are two successive Fibonacci numbers, hence L/L' = (1 + v5)/2 = T. In 
agreement with expected values shown in Tabled ui r = 1 (thus © r (av) = 1), and ui b = 3 (thus Q b (uj b ) — 3). For this 
geometry, rij = = 2. Plugging in these values into Eg. 1161 



J = {2-T)k 



1 3 

2 7+ 2^ 



pLL' « 0.2570 • pLL'. 



(17) 
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FIG. 10: Typical jam configurations. For type I jams, shown in (a), c = d — Tk' . For type II jams, shown in (b), d = d ma x and 
c = c max . Knowing the slopes of the lines in (a), s = l,s r = 2, and s& = 1/2, we can determine all of the angles: <j> = tan" (1/2), 
6 = tt/4, a = tan" 1 2, and £ = (tt/2 - tan" 1 2). 

Solving for the velocity, noting that the overall number of particles, N = pLL', 

v rl = 1 - — ps 0.7430. (18) 
N 

Note, this is independent of p and independent of L and L' . This predicted value for v r \ is included in the plot of the 
experimentally determined velocities, shown in Fig. [3] Note the calculated value slightly underestimates the number 
of particles involved in the jam (hence slightly overestimates v r ±). 

2. Rectangular aspect ratios, type II: 

As mentioned, rectangular lattices are well described by the formalism in Sec. IIVI We observe the "rich" jam 
described above, but also a second type of "depleted" jam (not enough particles). Empirically, all observations of this 
type have s ps 2/3, and one large jam. See for instance Fig. JSfc). The per row and per column width of the jams 
are the maximum, d max and c max respectively, since Tk! is greater than the number of particles available in a row or 
column. A jam interface of length k' and slope s (with s > 1) involves k! columns, but only sk' rows. Hence the total 
number of particles involved in the jams: 

J k (s • dxxxax C max ) 

= [(2u) r L' - L)/(2s-l)][spL/2u r e r (u r )+pL'O b (uj b )/2u b ], (19) 

For the realizations contributing to the plot in Fig. L/L 1 = (1 + y/S)/2 = T, io r — 1 (thus Q r (uj r ) — 1), ujb = 2 
(thus Ob(iOb) = 2), and s = 2/3 (which is empirically determined). Plugging these values in we find: 

J= l(j~ 1 )' + ' pLL ' ra °' 6293 ' pLL '- (20) 

Thus the average velocity, 

v r2 = 1 - — » 0.3707. (21) 
We include this predicted value in the plot of Fig. |SJ Note the agreement with the experimental data. 

3. Square aspect ratios: 

The typical geometry of an intermediate state for a square lattice is shown in Fig. iffii) . The interfaces and especially 
the edges of the bands are disordered and jagged. The slope assumptions only hold approximately: s r ps 1/2 and 
si, ps 2. Furthermore, there are several particles moving freely in the low density regions, unlike for the rectangular 
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lattices, where all particles eventually order into the bands and jams. Also unlike for the rectangular lattices (see for 
instance Fig.JSb)), the red and blue bands cross through each other without the pronounced shifting upwards. 

We cannot use the formalism developed above in Sec. lIVBl for this situation, since that formalism is based solely 
on geometric constraints of winding seamlessly around the lattice. A configuration on a square lattice with s r = 1/2, 
Sb = 2, s = 1, would have w r — 1/2. Plugging into Eqn. JfjJ) we find the required length of the overall interface, k' = 0. 
Instead, empirically we find the blue jams form a trapezoidal shape of approximate length b and height a/4. Likewise 
the red jams form a trapezoid of approximate length a and width 6/4. Each jam has this shape and there are rij = 3 
jams altogether (i.e., three distinct intersections of the bands). The number of particles involved in jams, J is: 

J « nj (ab/4 + ab/4) = 3a 2 /2. (22) 

Using Eq.^Jto solve for a and the fact that the overall number of particles, N = pLL', we can solve for the fraction 
of particles in the jammed state: 

J 3/3 r V 1 27 

Hence the velocity, 



N~2 U pL ) ^ = 32 P - (23) 



J 27 

v s ra 1 - 4r = 1 - —p. (24) 
N 32 p v ; 

This predicted value v s is included in the plots of Fig. [5J It captures the features of the experimental data, including 
the slight decrease in v s with increasing p. 



V. DISCUSSION 



The BML traffic model is a simple model of a jamming transition with self organization. In our study, instead 
of agreement with conventional beliefs, we find stable intermediate configurations with phase coexistence of jammed 
and free-flowing traffic. Such configurations have not been previously reported in the literature, despite the extensive 
amount of past work on the BML model. Furthermore, these intermediate configurations have interesting geometric 
and topological properties, with different behaviors resulting as a consequence of different aspect ratios of the underly- 
ing lattice. We develop a formalism, based on geometric constraints imposed by the lattice, to predict the asymptotic 
velocities of the coexisting phases. Visualizing the kinetic pathways of the evolving configurations was a key clement 
in uncovering the existence of the intermediate phases and, moreover, their periodic nature on lattices with relatively 
prime aspect ratios. The observations described in this manuscript open up a range of new questions about the BML 
model. 

As mentioned, instead of a phase transition as a function of density, we observe a bifurcation point where the 
intermediate states first begin appearing, and a second bifurcation point, where they completely cease to appear. 
Perhaps more interesting than predicting the asymptotic velocities, would be to calculate the locations of the bifur- 
cation points. From our experimental data, the exact location of the bifurcation points are difficult to determine, and 
moreover, also depend on the aspect ratio of the underlying lattice. 

It is possible that there is a sharp phase transition. However, in such a case, the density p would not be the 
appropriate order parameter. Perhaps a more appropriate order parameter would be an interaction energy between 
north-bound and east-bound particles. Note that when in the free-flowing state, the north and east particles have 
moved onto non-interacting lattices. It may be possible that one can define an initial energy based on the overlap or 
interaction between two lattices, and use that as an order parameter. 

A complication which makes theoretical treatment of the BML model difficult, is that it is not strictly monotonic. 
Adding particles to a configuration that is known to jam (i.e., increasing p), can actually change the sequence of 
particle interactions and result in that configuration going to free-flowing instead of jamming. Furthermore, it is 
known that certain discrete models with the same property as BML — namely that the randomness is in the initial 
condition, yet the dynamics fully deterministic — can be notoriously difficult to deal with analytically. Examples 
include bootstrap percolation[2l| and the Lorentz lattice gasj^. We modified the BML model to include a small 
probability for particles to flip species-types at each update. Our preliminary studies, adding this small amount of 
randomness to the dynamics, suggest that the model with randomness has extremely different geometric properties 
than the original BML model. In addition, for the model with randomness, we did not observe the intermediate 
configurations described herein. 
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